#include "GenThreads.h"
#include "loader.h"
#include "read_afof.h"
#include "Spline.h"
#include "Histogram.h"
#include "globvars.h"
using std::pair;
void  GetOneRot(vector<double> Time, std::pair<double, double> BarRot,
		double Tcurrent,double *Ts, double *Te, 
		vector<double> Phibar)
{

  //

  int Np=Time.size();
  int Np2=Phibar.size();
  CSpline<double> splineTtoPhi(
			 (double*)&Time[0],
			 (double*)&Phibar[0],
			 Np,"linear ");
  
 double Phicurr=splineTtoPhi.y(Tcurrent);

 double PhiStart=Phicurr+BarRot.first;
 double PhiEnd=Phicurr+BarRot.second;

 CSpline<double> splinePhitoT( (double*)&Phibar[0],
		     (double*)&Time[0],
		     Np,"linear ");
 *Ts=splinePhitoT.y(PhiStart);
 *Te=splinePhitoT.y(PhiEnd);
}

void CApVec::save_ap()
{
 
  ofstream fo(m_filename.c_str());
  
  for(tItApPeaksMap it=ApVec.begin();
	     it!=ApVec.end(); it++)
    {
      
      fo<<setprecision(4)<<std::fixed<<
	(*it).first<<setw(4)<<"  "<<
	       (*it).second.size()<<setw(10);
      tItApPeaks itP=(*it).second.begin();
      while(itP!=(*it).second.end())
	{
	  fo<<" "<<
	      (*itP).pos[0]<<setw(10)<<" "<<
	      (*itP).pos[1]<<setw(10)<<" "<<
	      (*itP).peak<<setw(1)<<" ";
	  itP++;
	}
      fo<<endl;
    }
  fo.close();
 
  // if(!FileExists(fraw))return;
  int np=Ap.size();
  
  ofstream of(fraw.c_str(), std::ios::out | std::ios::binary);
  of.write((char*)&np, sizeof(int));
  of.write((char*)&Ap[0], sizeof(float)*np);
  of.close();
  
}

bool CApVec::load_ap()
{
  return false;
	if(!FileExists(m_filename))
		{
		cout<<"Cannot find File with ApVec information"<<endl;
		cout<<"Failed to open: "<<m_filename<<endl;
		return false;
		};
	cout<<"# Reading Apocenter and pericenter peaks information filen: "<<m_filename<<endl;
	ifstream ifile(m_filename.c_str());

	ApVec.clear();
	unsigned int line=0;
	while(!ifile.eof())
		{
		 float x, y;
		 bool peakflag;
		 int snap, nr;// nr is the number of peaks in the given snapshot
		 ifile>>snap>>nr;
		 if(nr < 3)
		   return false;
		 for(int i=0;i<nr;i++)
		   {
		     ifile>>x>>y>>peakflag;
		     tPeaks val(x, y, peakflag);
		     ApVec[snap].push_back(val);
		   };
		 if(ifile.fail() && !ifile.eof()){
		   ifile.close(); 
		   cout<<"# some bad data in  ApPeaks: FILE LINE: "<<line<<endl;
		   return false;		   
		 }
		 line++;
		}	
	ifile.close();
	
	
	if(!FileExists(fraw))return true;
	ifstream irfile(fraw.c_str(), std::ios::in | std::ios::binary);
	int np;
	irfile.read((char*)&np, sizeof(int));
	Ap.resize(np);
	irfile.read((char*)&Ap[0], sizeof(float)*np);
	irfile.close();
	cout<<"# Reading done."<<endl;
	return true;
};

  bool   CApVec::GetAp(const int isnap, tApPeaks &c)
    {
      if(ApVec.find(isnap)!=ApVec.end())
	{
	  c.clear();
	  c=ApVec[isnap];
	  return true;
	}
      m_updated=true;     
      return false;

    };
void  CApVec::ApCom(const int isnap, const tApPeaks c)
  {
    
    ApVec[isnap].clear();
    ApVec[isnap]=c;
    m_updated=true;  
  };
/////////////////////////////

bool CComVec::load_com()
	{
	 
	if(!FileExists(m_filename))
		{
		cout<<"Cannot find File with bar COM information"<<endl;
		cout<<"Failed to open: "<<m_filename<<endl;
		return false;
		};
	cout<<"# Reading COM filen: "<<m_filename<<endl;
	ifstream ifile(m_filename.c_str());
	ComVec.clear();
	unsigned int line=0;
	while(!ifile.eof())
		{
		 float a, b, c, d, e, f;	
		 int snap;
		 ifile>>snap>>a>>b>>c>>d>>e>>f;		
		 if(ifile.fail() && !ifile.eof()){
		   ifile.close(); 
		   cout<<"# some bad data in  COM: FILE LINE: "<<line<<endl;
		   return false;		   
		 }
		 ComVec[snap]=CComPos(a, b, c, d, e, f);		 		 
		 line++;
		}	
	ifile.close();
	cout<<"# Reading done."<<endl;
	return true;
	}

void CGetData::GetComByMST(float *COM)
{
  string Warn="# we do not support this fueture.\n "
  "You should have look in to boostmst project";
  cout<<Warn<<endl;
  return;
}

template <class T>
struct setRhoGt {
  T& that;
  setRhoGt(T& f):that(f){}
  bool operator( )(const unsigned int p1,const unsigned int p2) {
    return(that.pRHO[p1] > that.pRHO[p2]);
  }
};
void CGetData::GetComByAFOF(float *COM, string afofmask)
{
 
  {
    //boost::mutex::scoped_lock scoped_lock(io_mutex);      
 
  CFOFCatalog *cat;

  char buf[512];
  sprintf(buf,afofmask.c_str(), m_isnap);

  string name=string(buf)+string(".grp.DAT");
  string name2=string(buf)+string(".grp");
  if(!FileExists(name)||!FileExists(name2))
    {
      GetCom(COM);
      return;
    }

  
  cat=new CFOFCatalog(name.c_str());
  cat->InitIDs( pID);
  cat->SetWeightForIDs( pRHO);
  double lCOM[]={0,0,0,0,0,0};
  unsigned int n=cat->GetNhalo(),i=0,idx=0;
  

  

  vector<unsigned int> setIndex; 

  if(n>0)
    {
      CGalaxy *gal=cat->GetGalByID(0);  
      n=gal->idref.size();
      setIndex.resize(n);      
      //      ofstream of("tst.txt");
      
      for(idx=0;idx<n;idx++)
	  {
	    setIndex[idx]=idx;  
	  }
      std::sort( setIndex.begin(), 
		 setIndex.end(),
		 setRhoGt<CGetData>(*this)
		 );
	
      unsigned int counter=0;
      for(vector<unsigned int>::iterator it=setIndex.begin();it!=setIndex.end();++it)
	{
	  i=gal->idref[(*it)]*3;

	  lCOM[0]+=pPOS[i];
	  lCOM[1]+=pPOS[i+1];
	  lCOM[2]+=pPOS[i+2];
	  lCOM[3]+=pVEL[i];
	  lCOM[4]+=pVEL[i+1];
	  lCOM[5]+=pVEL[i+2];   
	  //	  if(counter++ > n*0.9)break;
	}

      //of.close();

      	 
      for(i = 0;i<6;i++)
	{
	COM[i]=float(lCOM[i]*1.0/n);	
	}
      //      cout<<"SNAP="<<m_isnap<<" "<<COM[0]<<endl;
      // exit(0);
     
    }
  
  delete cat;
  }  

}
